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Introduction 

Calcium signaling is an intermediate step in many of the signaling pathways in neurons 
and glial cells and is informative of functional neural activity. In neurons calcium sig- 
naling precedes sub threshold and threshold (i.e. action potential) changes in membrane 
voltage, and can be used to infer electrophysiology from optical imaging ^21 17[ EHl l35] . 
In astrocyte glial cells it underlies the mechanisms by which these cells communicate in 
astrocyte networks and in bi-directional communication with neurons [l]|6l|30]. Relative 
changes in cytosolic calcium concentration can be measured using different fluorescence 
indicator dyes that can be imaged by optical microscopy in the visual light range, such as 
bulk loaded AM esters and genetically encoded calcium indicators j^SlSl]. The emitted 
fluorescence of indicator dyes change as a function of the relative amount of free calcium 
ions individual indicator molecules are able to interact with. Although the relationship be- 
tween measured fluorescence signals and the calcium levels that produce them is complex 
and non-linear, it is assumed that there exists a correlation between measured changes in 
emitted fluorescence by indicator molecules and differing cytoslic calcium concentrations. 
In this context, the measured fluorescence signal provides a valuable qualitative metric of 
changing calcium levels that allow inferences of cell signaling and function. Throughout 
the rest of this paper, we will use the terms "calcium signal" or "calcium fluorescence" 
to mean a measured calcium indicator fluorescence signal that reflects a relative cytosolic 
calcium concentration, as is routinely implied in the literature, even though in practicality 
we never know the real, i.e. absolute, free ion concentration that gives rise to the measured 
fluorescence signal. 

The data collected by a typical experiment records qualitative movies of imaged changes 
in calcium fluorescence intensity. One can visualize calcium transients and their relative 
positions and durations, but there is no inherent quantitative analysis of the data by the 
experiment itself that allows one to derive the dynamics that characterize such signaling 
events. For example, things such as propagation speeds and directions (i.e. velocity), 
the kinetics of measured waveforms, or analysis that depend on such properties, such as 
identifying and mapping the signaling geometry of intercellular calcium waves in networks 
of neurons or astrocytes. Measuring and tracing calcium (or other second messenger) 
fluorescence signals quantitatively from recorded movies manually is a tedious and labor 
intensive process for even small data sets, and involves comparing intensities at different 
frames and locations in order to calculate speeds and directions. It is generally not possible 
to do so for large data sets that encompass high spatial and temporal resolution detail or 
large numbers of cells interacting in a circuit or network. 

This can be addressed by analyzing experimental data with a filter algorithm called 
optical flow, which can be used to derive quantitative measurements of observed spatiotem- 
poral calcium signals imaged from fluorescence movies. The resultant vector data has a 
variety of uses, ranging from deriving basic measurements of signal velocity and direction, 
to characterizing and classifying spatiotemporal calcium dynamics between different ex- 
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perimental conditions. Optical flow is an imaging technique (i.e. a filter) that calculates 
a two-dimensional displacement field between two subsequent frames in a movie, based on 
the local spatial and temporal gradients of the two images. The optical flow filter origi- 
nated in the computer vision field, where it was designed to approximate object motion in 
time-ordered image sequences for applications like stereo disparity measurements, motion 
estimation, movie encoding and compression, and object segmentation [15]. The algorithm 
uses a computed local spatial and temporal gradient to approximate a displacement or 
flow vector at each pixel in the image. In both neurons and glial cells cytosolic calcium 
concentration changes manifest themselves as transient responses with a rapid increase, i.e. 
rising phase, followed by a kinetically slower decaying phase. This is because free calcium 
is cytotoxic and therefore kept at nanomolar concentrations in the cytoplasm under normal 
conditions. It is only transiently elevated followed quickly by its re-uptake or extrusion. 
Temporal changes are typically coupled to spatial changes as a signal propagates through a 
cell. Measured fluorescence changes then trace specific paths during periods of observation 
(c./. Fig. [T]). Calcium transients start at a particular location, travel in some direction 
at a specific speed and terminate at a different location. The typical kinetics of calcium 
transients in neural cells are particularly well suited to the computational requirements of 
the optical flow algorithm. 

We have successfully applied optical flow to calcium fluorescence movies and obtained 
displacement vectors that track the spatiotemporal progression of calcium signals. The 
filter works for calcium fluorescence data because calcium signals exhibit both spatial 
and temporal gradients. The computed vectors provide point estimates of the speed and 
direction of signals. Optical flow is ultimately an imaging filter that works on whole 
movies, much like edge filters and image segmentation filters are used in static microscopy 
[T2t [25t [13], and provides a novel and automated way of analyzing the spatiotemporal 
dynamics of calcium intracellular signaling in neurons and astrocytes. 

We begin by briefly reviewing the mathematics of the optical flow algorithm, describe 
how to solve for the displacement vectors, and how to measure their reliability. We then 
compare computed flow vectors with manually estimated vectors for the progression of a 
calcium signal recorded from representative astrocyte cultures. Finally, we applied the al- 
gorithm to preparations of primary astrocytes and hippocampal neurons and to the rMC-1 
Muller glial cell line in order to illustrate the capability of the algorithm for capturing 
different types of spatiotemporal calcium activity. We discuss the imaging requirements, 
parameter selection and threshold selection for reliable measurements, and offer perspec- 
tives on uses of the vector data. 

Optical Flow Algorithm and Computation 

In this section we briefly introduce the concepts and mathematics of optical flow, focusing 
in particular on our own implementation of the algorithm to the experimental data that 
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follows in the Results section. The theory behind the algorithm is well established and the 
interested reader is referred to a number of excellent texts on the subject (see for example 
[miinj). Optical flow is an algorithm that operates at the pixel level and calculates local 
displacement or velocity between time ordered image pairs. Optical flow (or equivalently 
image flow) is the perceived motion of an object in a field of view (e.g. by the human eye or a 
camera), defined as the "flow" or change in space and time of gray values at the image plane. 
It is an estimation of the motion field, which is the actual motion of the object in three 
dimensional space projected onto the image plane (i.e. what we would like to know). As 
long as the frequency of successive frames in an image sequence is shorter than the motion 
or displacement of the object of interest (in order to avoid confounding ambiguities in 
detecting the components of the motion caused by aperture and correspondence problems- 
see \T5*, T6]), the optical flow algorithm is able to track the motion of objects in the field 
of view as a function of changing gray scale levels, subject to appropriate constraints and 
minimizations. In other words, the algorithm assumes that any changes in gray values are 
due to the object moving, and that the irradiance of the object is constant from frame 
to frame. (This is actually a weak assumption that is difficult to satisfy since motion 
usually causes changes in irradiance, which is why the algorithm is an estimation of the 
motion field. In cases where irradiance does not change, the optical flow exactly equals the 
motion field.) The algorithm assumes the conservation of gray levels in the field of view 
and assumes that any changes in the distribution of gray levels are due to motion. In fact, 
the optical flow constraint equation (introduced below) can be derived by analogy from 
the continuity equation in fluid dynamics that conserves mass ([16j). By computing the 
optical flow for all pixels in a field of view, displacement vectors can be calculated for each 
pixel that map where an object moved to from the pixel in the first frame to that in the 
second. Intuitively, one can see why the algorithm performs best with objects that have 
strong contrasts at boundaries or large signal to background noise ratios. The kinetics of 
calcium transient signals display clearly distinguishable rising and decaying phases that 
trace specific paths during periods of observation in the form of intracellular calcium waves 
(Fig. [T]) that are readily detectable by the algorithm. 

The underlying assumption for computation is to constrain local temporal gradients 
to the product of spatial gradients and displacement vectors. The basic principle of the 
algorithm takes as inputs two images and computes a vector for each corresponding pixel 
in the images which approximates the displacement of a small window surrounding that 
pixel between the two images (Fig. [2]). Only intensity values inside the window are used for 
computing the pixel displacement value, so the measurement is localized. Adjacent pixels 
will have overlapping windows, so their computed vectors will be similar, much like pixels 
in a blurred image are similar. Following a mathematical description of the algorithm we 
describe the method for its solution and implementation that we used to derive the optical 
flow for calcium signals. We also discuss parameters and constraints of relevance to calcium 
fluorescence movies. 

Consider an arbitrary pixel with gray level intensity /(x, y, t), displaced in the xy plane 
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Figure 1: Selected frames from recorded movies of imaged calcium fluorescence activity in 
sparse networks of primary dissociated cortical astrocytes (top) and hippocampal neurons 
(bottom- where a single neuron at high magniflcation is shown). The color coded scale 
bars on the right represent fluorescence intensity / in units of Al/sec, as a flrst derivative 
of the calcium signal. For neurons the signaling was spontaneous, while for the astrocytes 
waves were pharmacologically induced. Fluorescence increases followed a relatively smooth 
spatial progression across the frames at the times indicated by the time stamp in the upper 
left hand corner of each image. Areas of increasing calcium concentration appear as positive 
Al/sec values, while areas of decreasing calcium concentration appear as negative values, 
but at a much smaller magnitudes. 
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Figure 2: The optical flow algorithm. A window Q in the same location in two subsequent 
image frames is used to compute a displacement or flow vector (arrow) for the pixel at 
the center of the window. Only image intensity values in are used for the calculation. 
Vectors are computed for each pixel in an image frame except in border regions where ft 
falls outside of the image. Given the position of the pixel as (x, y) at t seconds, (x, t), the 
displacement vector defines the motion of the pixel at the subsequent frame at 6t seconds 
as (x + 5x, y Sy^t -\- St). 
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by 6x and 6y at time 5t in an n x n window Q (Fig. [2]). This implies that 

/(x, y, t) = I{x + fe, y + 5y, t + 5t) (1) 

A first order Taylor series approximation of I{x,y,t) by expansion of the right side of 
[T] results in 

dl dl dl 

I(x + 6x, y + 6y, t + 5t) — I(x^ y, t) + — 5x + ^5y + —5t + higher order terms (2) 

ox oy ot 

Ignoring higher order terms, which provide negligible contributions, and taking into 
consideration equation [l] 

Dividing by 5t 

dl 5x dl 5y dl ^ , ^ . 

dl dl 91 ^ 

The two spatial and one temporal gradients are defined by and respec- 

tively. Ux — ^ and "^y — % represent the x and y spatial components of the optical flow 
displacement vector \x{x^y) — {ux,Vy). The basic optical flow formulation is to constrain 
temporal intensity changes (gradients) to the product of spatial gradients and u(x, y) to 
give equation [4j In more compact notation this can be written as 

VI{x,y,t).u{x,y) + ^^^^^^0 (5) 

Computing optical flow means finding the values of u(x, y) at each location for every time 
point that satisfy the above constraint, given the known local image intensity spatial and 
temporal gradients. 

Two factors establish comput ability of meaningful non-zero flow vector values. First, 
local spatial gradients must be non-zero at the point of interest {x,y,t). There has to 
be some image information around the pixel of interest, meaning that neighboring points 
have to have different values so that gradients are non-zero. If all pixels in a window 
around {x,y,t) have the same intensity values, then spatial gradients are zero and motion 
is undetectable by any means. Second, for displacement between subsequent frames to 
be computed, there has to be a temporal gradient at (x,?/,t), or some change in intensity 
between time points. If there is no temporal change in intensity between subsequent time 
points, then a value of u(x, y) = satisfies the constraint equation in [5j Both of these 
requirements are limitations on the original application of the optical flow when estimating 
displacement in natural scenes: objects may have constant intensity in a small window and 
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still be moving, meaning that motion may occur and the recorded intensity spatial and 
temporal gradients equal zero. These limitations are less important when optical flow is 
applied to calcium fluorescence movies. 

There are many methods for calculating optical flow for recorded movies (see [5j for a 
review), and all of them work on digitized movies with discrete pixel values of position and 
time, i.e. G {columns, rows^ frames). We chose the Lucas and Kanade method 

because it is conceptually simple and efficient, and flexible in terms of the image processing 
steps required for computation [3", "3^/22]. First, computation of the flow vector u{x,y) is 
performed on a window or spatial neighborhood of arbitrary size, centered around (x, y), 
which is more reliable than a single point estimate at Second, a window function 

W{x, y) is defined to favor values in the center over those near the edges. The choice 
of window size will depend on a variety of factors. It must be large enough to capture 
the apparent displacement across frames and small enough to resolve features of interest. 
The capture frame rate must be fast enough for displacements to be observable within the 
width of the spatial observation window across successive frames. When measuring the 
spatiotemporal motion of calcium signals the size of the window fi, the frame rate, and the 
resolution are all deeply tied to the size of the cells or cellular compartments in which the 
signal travels. Together, these parameters must be chosen so that the signal is observable 
and smooth enough to measure reliably as flow vectors across frames. For example, the 
choice of parameter values used to image calcium signals in part of a dendrite or a fine 
astrocyte process will necessarily be different than parameter values for broad calcium 
signals that fill the soma. The constraint equation is redefined as a weighted least-squares 
fit of local first-order constraints to a constant model of a local u(x, y) in each small spatial 
neighborhood around the pixel of interest. The goal is to find the value of u{x,y) that 
minimizes 



J2 W''{x,y)(vi{x,y,t)-u{x,y) + 



dt 



The above equation can be rewritten and solved as the linear system: 

A^W'^A ■ u(x, y) = A^W^h 



(6) 



(7) 



Where, for neighborhood f2, consisting of n points centered around the pixel and time of 
interest ix,y,t), Q = {{xi,y2,t), {x2,y2,t), {xn,yn,t)}: 



A 



%ixi,yi,t) ^{X2,y2,t) ... %ixn,yn,t) 

91 +\ dl(^ 4.\ dl I 



d^{xi,yi,t) Q^{x2,y2,t) ... Q^{Xn,yn,t) 

W = diag\W{xi,yi), . . .,W{xn,yn)] 



b = - 



dl dl 1 ^ 

— . . . , -^{Xn,yn,t) 



(8) 

(9) 
(10) 
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ft is usually a square window with sizes typically ranging from 3x3tol5xl5orn = 9 
to n = 225 points. We have set the weight matrix to a two dimensional Gaussian with 
equal to 1/6 of the window width. As an example, for a 5 x 5 or n = 25 point window: 
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Here, the center values in W have a greater contribution to the calculation than the edge 
values, favoring gradient values at the pixel of interest. 
Solving for the flow vector u{x,y) in equation [?[ yields: 



u{x, y) = [^'^1^2^] ^A^W^h 



(11) 



Equation [TT] describes a linear system in matrix form, where the flow vector u at spatial and 
time location (x, y, t) is solved from the quantities of A, VF, and b, defined from the spatial 
and temporal derivates of n points around [x^y^t). The 2x2 matrix [A-^VF^A] matrix 
contains all the image spatial derivatives, and if those values are close to zero, the matrix is 
poorly conditioned, and flow estimates become unreliable. Ensuring that both eigenvalues 
of the [A-^VF^A] matrix are sufficiently large is a good way to ensure that the matrix is 
well conditioned, since a measure of the conditioning number is the ratio of the largest to 
the smallest eigenvalue [5j. While this is not the only way to compute conditioning, this 
is the test we used for visualization and measurement reliability of computed vectors for 



calcium fluorescence data (see Appendix A. 3 for more information 



Spatial and temporal derivates were computed using 2x2 convolution kernel filters, 
where the ** operator denotes 2-dimensional discrete convolution: 



dl{x,y,t) 1 
dl{x,y,t) 1 

dl{x,y,t) _ l.(j(^^y^i + ^i)_j(^^y^i)),,l 



-1 1 

-1 1 

-1 -1 

1 1 



dt 



At 



1 1 
1 1 



(12) 
(13) 
(14) 



Here At represents the time between frames or the frame rate 1/At. Since the temporal 



derivative calculated in [M] forms the basis for the b vector in[TT} the frame rate has a linear 
effect on the magnitude of the flow vector u. 

Optical flow outputs a displacement vector in units of pixels, normalized to the time 
difference between the two frames used for computation. When normalized, the vector 
takes on velocity units of pixels per frame (for this reason it is called a flow vector). 
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The conversion to physical units will depend on the spatial resolution of the camera and 
microscope, typically expressed in microns per pixel, and the sampling rate for the movie 
capture, expressed in frames per second. Spatial resolution is a function of the objectives 
used as well as the resolution of the imager and any pixel binning used. The frame rate is 
limited at the high end by the camera sampling rate, and at the low end by the minimum 
exposure time required to capture a detectable intensity signal. The exposure time may be 
reduced by increased gain or pixel binning, but those come at a cost of reduced resolution or 
increased noise. The conversion between units of pixels/frame and units of microns/second 
is straightforward: 

microns pixels frames microns 
second frame second pixel 

While the optical flow algorithm produces vectors in units of pixels/frame, the analysis 
of the data in the Results section below have been converted into physical units of mi- 
crons/second, using the resolution and capture rate of the recordings given the specifics of 
our imaging system. 



Results 

Comparison between computed and manually estimated flow vectors 

We manually estimated flow vectors for 12 images equivalent to 6 seconds of calcium sig- 
naling in primary dissociated spinal cord astrocyte cultures (orange arrows in Fig. [s]), and 
qualitatively compared them to computed optical flow vectors for the same data (green 
arrows in Fig. |3j note that only reliable vectors are shown as determined by the eigenvalue 
test- c.f. equation pT| and above discussion). Manual estimation required stepping through 
frames and approximating roughly how a calcium signal progressed in time, which in this 
experimental preparation included intercellular calcium waves that propagated through a 
subset of the cell network. The manually traced signals were not the only ones observ- 
able in the small movie sequence used, but were chosen to illustrate four representative 
signaling paths. Manual estimation was performed in two second intervals, estimating the 
incremental spatial progression of a given calcium signal across four frames. 

Estimation of the flow or displacement of a cell signal such as calcium between frames 
manually like we did for the data in Fig. |3] is a very tedious and labor intensive process, 
and can only realistically be done under very sparse conditions where the observer can 
clearly delineate the flow of the signal visually. It is nearly impossible to do at the pixel 
or small window level. In contrast, optical flow calculates a displacement vector for every 
pixel in every frame, operating at a much finer scale and capturing much more detail than 
is possible with manual estimates. Nonetheless, in Fig. |3]for the purpose of qualitatively 
validating derived optical flow vectors to manually estimated ones, in both cases there was 
a clear overlap in vector direction between manual and flow vectors. By contrast, there 
were greater differences in vector magnitudes between manual and flow vectors, which is 
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Figure 3: Comparison between computed optical flow vectors (green) and manually es- 
timated flow vectors (orange) for a primary dissociated spontaneously forming astrocyte 
network in culture. While computed vectors were calculated for every pixel and every 
frame, manual vectors were estimated every four frames and only trace a few selected sig- 
nals. Only reliable optical flow vectors are shown, and only one in four vectors in both 
horizontal and vertical directions are shown for clarity. Unlike the manual vectors, flow 
vectors are only shown for the current frame. The top sequence of panels show vectors 
overlaid on extracted frames from the actual movie at the indicated times for the entire 
field of view. The bottom panels show the vectors in detail for the 2, 4, and 6 second frames 
in order to more clearly assess the qualitative overlap between optical flow computed and 
manually estimated results. For optical flow vectors (in green) only vectors that putatively 
correspond to manual vectors (in orange) are shown, in contrast to the upper panels which 
show all computed vectors (see text). See appendix for details regarding experimental 
preparations, imaging, and parameters for calculation. 



11 



Buibas et. al. 



Mapping neural cell networks with optical flow 



consistent with the fact that the manual estimates spanned four frames while optical flow 
vectors were calculated across adjacent frames. There is also temporal overlap between 
the two cases, in the sense that similar displacements were estimated for the same frames 
using both approaches. The eigenvalue threshold masked out unreliable vectors, and this 
correlated well with calcium activity; only areas of spatial and temporal changes in the 
movie produced reliable vectors as assessed visually, which is ultimately the most accurate 
estimator of complex motions, but only if given the right conditions (e.g. conditions that 
allow the human eye to separate motion). The optical flow algorithm however, is able to 
provide reliable quantitative measurements of signaling dynamics at spatial and temporal 
scales simply not measurable by qualitative visual inspection or manual estimations of the 
data. 

Optical flow characterization of interceUular signaling 

We applied the optical flow algorithm to typical calcium fluorescence movies of sponta- 
neously forming sparse networks of neural glial cells and neurons in culture, and looked 
at the dynamics of intercellular signaling following pharmacological or mechanical stimu- 
lation. We purposely chose sparse networks because it facilitates the visual interpretation 
of the entire resultant vector field, but the algorithm itself can operate on any data that 
displays an appropriate signal. We recorded movies from the rMC-1 Muller glial like cell 
line, which mechanistically displays calcium signaling similar to Muller retinal glial cells in 
vivo [39], primary dissociated spinal cord astrocytes, and primary dissociated hippocampal 
neurons. Intercellular calcium waves in rMC-1 cells and astrocytes were mechanically in- 
duced by gently poking an initial cell without penetrating the cell membrane, while calcium 
waves in neuronal networks were induced by the localized pharmacological application of 
glutamate to one or a small group of cells (see the appendix below for details about exper- 
imental preparations and imaging parameters). In particular, intercellular calcium waves 
in astrocytes and related anatomically specialized macroglial cells such as Muller cells in 
the neural retina or Bergman glia in the cerebellum have been known to occur under ex- 
perimental conditions for several years now, and have recently been shown in vivo under 
both physiological and pathophysiological conditions in different parts of the brain, medi- 
ated by intracellular calcium transients that induce paracrine signaling, primarily through 
adenosine triphosphate (ATP) [191 |20l |14] . Astrocyte and related macroglial cells engage 
in bi-directional chemical signaling with neurons and have the ability to modulate and 
directly participate in information processing in the brain, which necessitates more than 
just interactions between neurons and almost certainly involves astrocytes somehow. The 
functional roles of glial intercellular calcium waves and their contributions to modulating 
neuronal information are not yet known, and in fact the dynamics of these signaling events 
and the conditions under which they occur are just beginning to be explored. 

The key parameter for computing optical flow using the Lucas-Kanade method is the 
window size fi, specified as a square of a given width (see above). It defines the local 
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neighborhood of pixels along a point of interest that is used to compute the spatial and 
temporal gradients required for the calculation. Though not required for computation, a 
minimum value for the eigenvalues for the matrix A^W^A should be specified to mask 
out unreliable measurements. This ensures that only reliable displacement vectors are dis- 
played and used for analysis. Since the intensity values are a function of the experimental 
setup, microscope, and camera, the A^W^A matrix and its eigenvalues will scale accord- 
ingly. The selection of the eigenvalue threshold is thus arbitrary, much like the selections 
of the camera gain, exposure time, and other imaging parameters are made to generate 
easily visible intensity values (see Appendix A. 3 for more information on selecting suitable 
eigenvalue thresholds). Table [l] shows the window sizes, eigenvalue thresholds, and capture 
frame rates used to calculate the vector fields shown in Fig [Ij The displacement vectors 
can be converted into velocity by equation [T5j The original calcium fluorescence movies 
and Matlab code written to implement the optical flow algorithm are freely available by 
contacting the corresponding author. 



Table 1: Image capture and optical flow parameters for shown figures 



Parameter 


rMC-1 Cells 


Astrocytes 


Hippocampal 
Neurons 


Frame Capture 
Rate (Hz) 


16.4 


8 


4 


Window Size 

(pixels at 
1.3/xm/pixel) 


11 


9 


11 


Minimum 
Eigenvalue - 

(Ai,A2) 
greater than 


11 


1.4 


0.3 



Neuronal cultures displayed derived optical flow vectors along processes as the calcium 
signal propagated throughout the network. As expected, computed vectors and the resul- 
tant vector field followed the geometry of connected processes (i.e. axons and dendrites) 
in the sparse network (Fig. [4^). The pattern of activation in this example proceeded 
diagonally from the site of stimulation in the upper left hand corner of the field of view. 
Some neurons activated at considerably longer times following the stimulus (i.e. out to 
7 or 8 seconds) most likely due to recurrent feedback signaling in the network which can 
last several seconds. Note how since only reliable vectors are plotted, as determined by 
the eigenvalue test, there are spatial discontinuities in the temporal progression of mapped 
signals, which reflect areas where the algorithm could not compute reliable vectors given 
the measured data. This may be especially true at lower magnifications as in the example 
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shown here for comparatively large fields of view that capture many cells. This represents a 
challenging task for the algorithm. Nonetheless, both the spatial and temporal progression 
of calcium signals are easily visible. The computed data, being in vector form, can com- 
plement existing methods like cross-correlation that use only cell body data to establish 
relationships between cells for example. 

Signal flow patterns were also computed for astrocyte and rMC-1 glial networks (Fig. 

and c, respectively). Astrocyte signaling showed rapid burst like radial patterns that 
was mostly complete by 2 seconds, with some smaller regions of cells activating later as far 
out as 6-7 seconds. This is consistent with descriptions of intercellular calcium waves re- 
ported previously [HI [261 El |3T] . rMC-1 cells showed qualitatively similar radial patterns of 
activation, with signaling occurring within about 3 seconds following stimulation. However, 
unlike the astrocyte response, where there was uniform signaling across the network near 
the site of stimulation, rMC-1 cells showed more heterogeneity in spatial activation pat- 
terns, with distinct clusters of cells activating and spreading calcium waves. The distances 
traveled by the waves in the rMC-1 example roughly agree with previous quantitative char- 
acterizations of calcium waves in similar preparations, on average displaying wave distances 
of about 60 /im over the first two seconds or so and distances between 50-100 /xm over about 
4 seconds [39] . It is interesting however that the spatial progression of the calcium signals 
in this example was not linear as a function of time, in the sense that cells roughly equidis- 
tant from the site of stimulation activated at different times, within about one second for 
some versus as late as three seconds for others. The relationship and dynamics between 
the spatial versus temporal properties of such waves are difficult at best and usually not 
possible to determine by visual inspection of recorded movies alone, and are not captured 
by calculations such as the one dimensional signaling speed of a progressing wave front. 
Furthermore, speed and distance calculations of neuronal and glial signaling across net- 
works of cells are usually corse approximations computed using low magnification movies 
that provide a sufficiently large field of view. In contrast, optical flow provides reliable 
single pixel vectors for any sized region of interest that represent very fine grain detailed 
descriptions of calcium signal propagation difficult to achieve otherwise. For the astrocyte 
data from Fig. [4]3, Fig. [5] illustrates the distribution of signaling speeds in /im/second for 
65 optical flow vectors for a small 10 x 10 pixel region equivalent to a 13 x 13 fim region 
in the field of view (orange box in the figure). Any region size of interest anywhere in the 
imaged field that might be of functional interest to the investigator can be similarly char- 
acterized. By way of rough comparison, optical flow calculated speeds for calcium signals 
in computed window were distributed from 1-10 fim/sec, and are roughly similar to those 
reported previously using more approximate methods, in the range of 5 to 10 /xm/second 
[HI |26l 121 [31]. The bimodal distribution in the figure reflects what is visually apparent 
in the source movie: some of the areas in the orange pixel region exhibit spatiotemporal 
displacement while others do not, indicating that calcium concentration changes propagate 
in specific regions with specific patterns. Manual estimates from the literature typically 
look at maximum propagation speeds, as seen in the second peak at about 9 ^m/second. 
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Figure 4: Computed optical flow vectors for induced calcium signals in spontaneously 
forming in vitro networks of (a) primary hippocampal neurons, (b) primary spinal cord 
astrocytes, and (c) the rMC-1 Muller glial-like cell line. Six frames from each representative 
recorded movie are shown with the computed vector field superimposed at times indicated 
by the time stamps in each frame (left set of six panels). Right panels: Composite temporal 
projections of the entire movies. The vector fields show the full spatial progression for the 
evolving calcium signals, with time (i.e. temporal progression) color coded by the color 
map (in seconds). Plotting the vector fields in this way allows the full spatiotemporal 
propagation of derived signals from entire jigovies to be summarized in a single image. 
This facilitates the qualitative visualization and identification of complex dynamic signaling 
patterns that would be difficult to detect otherwise, such as for example by simply "playing 
back" the movie. 
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Figure 5: Optical flow velocity magnitude distributions for the astrocyte data from Fig. 
^p. The flow vector magnitudes for reliable measurements in a 10 x 10 pixel (13 x 13 /xm) 
region (orange square) are shown as histograms. 



Discussion 

We describe and show the application of optical flow gradient methods for identifying and 
spatiotemporally mapping functional calcium signaling in networks of neurons and glia. 
Although we focused on networks of cells here, the method can be equally applied to the 
analysis of spatially detailed sub-cellular compartmentalized regions of interest, such as 
dendrites or astrocyte processes. The method makes use of the spatial first derivative of 
moving objects in a field of view, in this case changes in fluorescence levels of calcium 
indicator dyes associated with the free concentration of intracellular calcium, to track their 
motion between subsequent frames in an image sequence (i.e. a recorded movie). The 
mathematical foundations of optical flow are well established and optical flow algorithms 
have been used in a wide variety of fields including applications to cell and molecular biol- 
ogy to track the movement of proteins, vesicles, and even whole cells [25 . In neuroscience 
and neural engineering it has been used in electromyography [181 123 sensory percep- 
tion [214 i28j, while clinically it has been used to detect seizures in neonatal infants [17j, 
among other applications. However, the method has not been previously applied to track- 
ing and visualizing calcium signaling and deriving quantitative measurements of calcium 
spatiotemporal changes that underlie intracellular and intercellular functional signaling in 
neural cells. 

Although in this paper we applied the optical flow algorithm to two dimensional flu- 
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orescence movies, the algorithm itself can be readily applied to a recorded movie made 
up of three dimensional stacks acquired using two-photon microscopy. Work by others is 
pushing two photon imaging towards recording real time functional signaling from three 
dimensional volumes of active cellular neural networks [TT|^j. If the sampling rate is suf- 
ficiently high, optical flow can be computed in three dimensions using a volume instead of 
a square window around a pixel to generate a three dimensional displacement vector. The 
same constraints on volume size, sampling, and vector reliability metrics in two dimensions 
apply to the three dimensional case. 

Optical flow methods produce a lot of data, generating a vector for every pixel in every 
image pair computed, so further processing, rendering and visualization methods are key to 
making quantitative comparisons between experimental setups. Statistical comparisons can 
be made from vector values by comparing differences between selected regions in different 
preparations; velocity averages for each region can be compared using statistical methods 
such as means, standard deviations, and p-values. While vector values from adjacent pixels 
are not statistically independent, averaged vector values for a given region of interest may 
be used for statistical comparison with another, non-overlapping region. 

Another potential use of the vectors is to classify spatiotemporal patterns. Similar to 
using a scalar kernel filter to match an image pattern such as an edge or corner, vector 
fields themselves can be filtered with a known vector kernel to match a pattern of interest. 
This method is called Clifford convolution [9J and has been used to label physical flow 
regimes in fluid dynamics applications. By designing a vector field filter and convolving 
it with computed optical flow vectors, a scalar map identifying specific patterns of flow 
associated with the saptiotemporal dynamics of the measured signal can be constructed in 
order to classify regions exhibiting such patterns. 

One of the most exciting potential uses of computed flow vectors is in functional network 
reconstruction. Borrowing again from the field of fluid mechanics, a dynamic vector field 
can be used to reconstruct the path of a hypothetical particle from a given starting point, 
tracing out the path that a signal might take between cells, much like a particle in a dynamic 
flow field [36j 37j . Geometrically mapped paths of measured signals that originate in an 
activating cell and propagate through a network may be very useful for reconstructing 
the dynamics of the network. This would complement existing network reconstruction 
algorithms which typically rely on temporal data around fixed regions of interest. 
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A Appendix 

A.l Cell Preparations 

rMC-1 glial cells and primary spinal cord astrocyte cultures (the latter dissected and grown 
similar to previously described [32^ |23j were grown to approximately 80% confluency and 
washed twice with Kreb-HEPES buffer (KHB) solution (10 mM HEPES, 4.2 mM NaHCOa, 
10 mM glucose, 1.18 mM MgS04, 7H20, 1.18 mM KH2PO4, 4.69 mM KCL, 118 mM NaCl, 
1.29 mM CaCl2, pH 7.4) and incubated with 5/xM Fluo-4AM in KHB for 1 hr at room 
temperature. Excess dye was removed by washing twice with KHB and an additional 
incubation of 30 min at room temperature was done to equilibrate intracellular dye con- 
centration and ensure complete intracellular hydrolysis. For astrocytes and rMC-1 cells, 
calcium transients were initiated by mechanical stimulation of a single cell using a (0.5/xm 
i.d.) micropipette tip (WPI Inc., Sarasota PL) mounted on a M325 Micrometer Slide Mi- 
cromanipulator (WPI Inc., Sarasota PL). Comparable data were obtained using adenosine 
triphosphate (ATP) pharmacological stimulation. 

Por hippocampal cultures, dissociated hippocampal neurons from timed-pregnant em- 
bryonic day 18 (E18) Sprague-Dawley rats were cultured on glass bottomed tissue culture 
dishes coated with poly-D-lysine and laminin (BD Biosciences, San Jose, CA). Cultures 
were plated at a cell density of 10^ cells/3. 8cm^. Cultures were maintained at 37C in 5% 
ambient CO2. Plating media was composed of basal medial Eagle (Invitrogen, Carlsbad, 
CA) with IX Glutamax, 1000 U/mL penicillin and streptomycin sulfate, 5% PBS, and IX 
N2 supplement. Culture media consisted of Neurobasal (Invitrogen, Carlsbad, CA) with 
IX Glutamax, 1000 U/mL penicillin and streptomycin sulfate, 20mM glucose, and IX B27 
supplement. Culture media was supplemented with lOuM Ara-C for 24 hrs at IDIV to 
inhibit overgrowth of glia. All imaging was performed on 3-5DIV. 

Bulk loading of hippocampal cell cultures was accomplished via incubation in the dark, 
at room temperature, for 30 min in 1/xM of the fluorescent Calcium indicator Pluo-4- 
AM in Krebs-HEPES buffer (lOmM HEPES, 4.2 mM NaHCOa, lOmM glucose, 1.18mM 
MgS04-7H20, 1.18mM KH2PO4, 4.69 mM KCl, 118mM NaCl, 1.29 mM CaCb, pH 7.4), 
followed by 2x 5 min washes in Krebs-HEPES with lOO/xM sulfinpyrazone. Hydrolysis 
was allowed to proceed for an additional 30 min. Stimulation of neurons with glucose was 
performed by microinjection of lOOuL of lOmM glutamate in PBS from a specified-side 
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of the culture dish, well outside of the microscope field of view. The fluorescence signal 
generated across the monolayer of cells was recorded for 10 sec prior to glutamate injection, 
and for 120 sec following injection. Cultured neurons were incubated for 30 min prior to 
imaging in Mg^+-free PBS to induce the synchronization of calcium transients. 



A. 2 Imaging Setup 

Visualization of calcium indicator dye fluorescence was achieved using a 488 nm (FITC) 
filter on an Olympus 1X81 inverted fluorescence confocal microscope (Olympus Optical, 
Tokyo, Japan) that included epifluoresence, confocal, phase, brightfield, and Hoffman dif- 
ferential interference contrast (DIG) modalities. Real-time movie recordings of calcium 
transient propagation were acquired with a Hamamatsu ORCA-ER digital camera (Hama- 
matsu Photonics K.K., Hamamatsu City, Japan) and Image-Pro Plus data acquisition and 
morphometric software (version 5.1.0.20, Media Cybernetics, Inc., Silver Spring, MD) or 
Lab View custom written data acquisition software (ScopeController). All images were 
captured with a lOX objective, using a 2x2 binning on the camera, for a resolution on 
1.3/xm/pixel and a total image size of 612x572 (camera's maximum resolution is 1224x1144). 
Images sampled at frequencies ranging from 2 to 16.4Hz, or 0.5sec to O.OGsec exposure time. 



A. 3 Reliable Vectors via the Eigenvalue Test 



Recall that flow vectors are computed from the linear system in 11 , This is a typical linear 
system of the form: 

M • u = z 

where u is the unknown and z and M are known quantities. The condition number of a 
matrix simply describes how a small deviation in the known z translates to an error in 
u. A high condition number means the matrix is ill-conditioned, meaning that a small 
deviation in z leads to a large deviation in u, making that computation unreliable. One 
way to compute the condition number of a matrix is to take the ratio of the largest to 
smallest eigenvalue of that matrix: 



Since M is a 2 x 2 matrix, it has 2 eigenvalues so ensuring that both are above a certain 
value makes the condition value relatively low. The minimum threshold value depends 
on the incoming intensity values. 

Intensity readings from the CCD camera can take on any number of values, based on 
the digitization (8-bit, 12-bit, 16-bit, for example), the exposure time, gain setting on the 
camera, and above all the dye loading in the cell preparation. Typically, during observation, 
the experimenter manually adjusts gain and exposure time to obtain reasonable intensity 
values, typically in the middle of the digitization range. 



22 



Buibas et. al. 



Mapping neural cell networks with optical flow 



Eigenvalues for the A^WA used in flow vector calculation typically scale with the range 
of recorded intensity values and are calculated for every pixel, producing an eigenvalue 
image map. The values chosen in table [l] were manually chosen during examination of the 
minimum eigenvalue image for a few representative frames, ensuring that they fell between 
ares where we visually detected spatiotemporal changes in intensity and areas were we did 
not detect such changes. This is the same process one would undertake when thresholding 
a regular monochrome image for the counting of cells: the intensity threshold is set to a 
value between the intensity of an area where there is a cell and an area where there is no 
cell. 
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